Predicting crystal structures: the Parrinello-Rahman method revisited 
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By suitably adapting a recent approach [A. Laio and M. Parrinello, PNAS, 99, 12562 (2002)] 
we develop a powerful molecular dynamics method for the study of pressure-induced structural 
transformations. We use the edges of the simulation cell as collective variables. In the space of 
these variables we define a metadynamics that drives the system away from the local minimum 
towards a new crystal structure. In contrast to the Parrinello-Rahman method our approach shows 
no hysteresis and crystal structure transformations can occur at the equilibrium pressure. We 
illustrate the power of the method by studying the pressure-induced diamond to simple hexagonal 
phase transition in a model of silicon. 

PACS numbers: 61.50.Ks, 64.70.Kb, 02.70.Ns, 07.05.Tp 
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Predicting equilibrium crystal structures at a given 
pressure and temperature is an important problem in 
fields of science as different as solid state physics, geo- 
physics, planetary physics, materials science, polymer 
science, chemistry, etc. Upon increasing external pres- 
sure crystals usually undergo structural phase transi- 
tions. Often, the final structure is unknown and sim- 
ulations can be very useful in identifying possible candi- 
dates. This task represents a challenge for computational 
physics. Great progress was achieved with the intro- 
duction of constant-pressure molecular dynamics (MD) 
|jj and in particular the Parrinello-Rahman method ||] 
where the box is allowed to change its shape in order to 
comply with a new structure. The Parrinello-Rahman 
method is now described in textbooks and widely used 
also in different variants § f§ | |, @, |, § @. How- 
ever, structural transitions are often first order. Since 
for crystal simulations periodic boundary conditions are 
commonly used, heterogeneous nucleation is suppressed 
and the system has to cross a significant barrier in or- 
der to transform from one structure to another. As a 
consequence large hysteretic effects are observed within 
the above approaches. In order to observe the transition 
within the accessible simulation time|ll[] one often has to 
overpressurize the system close to the point of mechanical 
instability. Under such conditions one or more interme- 
diate phases may be skipped^, [l2|, which reduces the 
predictive power of the method. 

In the work of Parrinello and Rahman|2| it was realized 
that in an MD simulation of a crystal phase transition it 
is necessary to treat the MD supercell edges a, b, c as 
dynamical variables. These variables were arranged to 
form a matrix h = [a, b, 5) and extending the work of 
Andersen ffl a Lagrangian was introduced that coupled 
the h degrees of freedom with the microscopic motion of 
the atoms under condition of constant pressure. Due to 
the time-scale problem mentioned earlier this approach 
tends to be ineffective at pressures close to the critical 
transition pressure. As the origin of the problem is the 



lack of efficiency of standard molecular dynamics in cross- 
ing high barriers, we adopt here a conceptually different 
strategy. Since our aim is to simulate a phase transition 
at a pressure P and a temperature T we consider the 
Gibbs potential Q(h) — J-(h) + PV as a function of h 
where J-(h) is the Helmholtz free energy of the system 
at fixed box and V = det(h) is the volume of the box. 
We assume, following Nose and Klein jj^], that the matrix 
h is symmetric in order to eliminate rotations of the su- 
percell. This reduces the number of collective variables 
to 6. These 6 independent components of h now repre- 
sent collective coordinates, or order parameters, which 
distinguish the different minima of Q. We note that the 
derivative 
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where p is the internal pressure tensor, can be easily eval- 
uated from microscopic MD or Monte Carlo runs at con- 
stant h by averaging the microscopic virial tensor. Mak- 
ing use of Eq.(j^) we now construct an algorithm, based on 
the recently introduced method of Ref . |p| , that is able 
to explore the surface £/(h) efficiently and in particular 
can identify the local minima which correspond to stable 
or metastable crystal structures at a given pressure P. 
The method jl3| has been shown to be able to dramati- 
cally speed up the simulation of activated processes and 
is therefore well suited for simulating first-order phase 
transitions. We describe here the basic ideas and refer 
for more details to the original paper. 

Following Ref.Jll|, the collective variables that are 
now arranged to form a 6-dimensional vector h = 
(^11,^22,^33,^12,^13,^23) are evolved according to a 
steepest-descent-like discrete evolution with a stepping 
parameter Sh (metadynamics) 
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Here, the driving force 0* = —75c- is derived from a 
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history-dependent Gibbs potential 5* where a Gaussian 
has been added to 5(h) at every point h* already visited 
in order to discourage it from being visited again. Hence 
we have 

t' 2 
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and the force <fr is therefore a sum of a thermodynamical 
driving force F — an d the term F g coming from a 
potential constructed as a superposition of Gaussians. As 
time proceeds the history-dependent term in Eq.([|) fills 
the initial well of the free-energy surface and the system 
is driven along the lowest free energy pathway out of the 
local minimum. The passage through the transition state 
can be detected by monitoring the relative orientation of 
the forces F and F g . While a well is being filled these two 
forces approximately balance each other, F + F g £3 0, and 
the two vectors have roughly opposite directions. After 
crossing the saddle point this is no longer true and F and 
F g become almost parallel and oriented along the eigen- 
vector corresponding to the negative eigenvalue of the 
Hessian matrix. The indicator F • F g /(|F||F g |) develops 
a sharp spike which can be used to signal the transition 
from one basin to the other. 

The choice of the parameters W and Sh depends on 
the 5(h) landscape. In order to achieve the necessary 
energy resolution, W should be chosen as a fraction of 
the relevant energy barriers. The parameter Sh on the 
other hand determines the resolution in h. However, a 
very small value of Sh is not to be recommended. In 
fact a small Sh requires longer runs. Furthermore for an 
optimal filling the curvature of the Gaussians should be 
smaller than that of the well. This leads to the condi- 
tion < x where K is the smallest eigenvalue of the 
5(h) Hessian at the minimum ho- For a cubic system we 
can estimate K from the approximate expansion of 5(h) 
around h 

5(Ah)«5(h ) + iyc(^) 2 (4) 

where L is the cell edge and c is of the order of magni- 
tude of the elastic constants. This leads to the estimate 
K ss Lc and to the condition -W? < Lc. A more general 
discussion of the choice of W and Sh can be found in 
Ref.||. 

In practice the metadynamics simulation proceeds as 
follows. We start from an equilibrated value of h at a 
given pressure P and temperature T and evaluate the 
pressure tensor p in a constant h MD run long enough 
to allow relaxation to equilibrium and sufficient averaging 
of p. The h is updated using the forces (|l|) and metady- 
namics equations (§,||) to a new value h . After the box 
is modified the particle positions are rescaled in order to 
fit into the new box using the relation r = h h _1 f. As 
the initial free energy well is gradually filled the box un- 
dergoes a set of deformations until a transition state is 



reached and the system enters into the basin of attraction 
of a new state. In order to characterize the new phase 
it is often useful at this stage to switch off the Gaussian 
term, so that the metadynamics becomes purely steepest 
descent-like and drives the system towards the equilib- 
rium state for the new structure. In this equilibrium 
state the pressure will be equal to P. However, dur- 
ing the metadynamics the pressure tensor can become 
anisotropic and the internal pressure may be different 
from P. Once the new structure is characterized one 
can switch the Gaussians on again, thus filling the new 
minimum, and move to other minima, if available. The 
metadynamics is capable of reconstructing the free en- 
ergy profile []l3]] , since the sum of the Gaussians in Eq.(||) 
converges to —5(h) up to an additive constant, if W and 
Sh are properly chosen. This will not be used here since 
once the structures are known it is relatively straight- 
forward to calculate their free energy |Q. We empha- 
size, however, that free-energy calculations alone do not 
provide an alternative to our method since they assume 
knowledge of the final crystal structure. 

We have tested our method on several model Hamilto- 
nians. Here we report only a simulation of a tight-binding 
model of SijljJ at a pressure very close to the theoretical 
transition pressure. This tight-binding parametrization 
captures some of the main features of the Si phase dia- 
gram and provides a convenient test model. In the follow- 
ing we shall use a supercell of 216 atoms and only the T 
point of the BZ. The T — phase diagram for this model 
system can be found by performing energy versus volume 
calculations in the three relevant structures, namely the 
P = equilibrium diamond structure and the two high- 
pressure phases, /3-tin and simple hexagonal (SH). The 
latter two are almost degenerate in energy, /3-tin being 
only metastable. A common tangent construction gives 
a critical pressure of 15.5 GPa for the transition from 
the diamond to the SH phase. Applying the Parrincllo- 
Rahman method to the same model and system size, the 
transition from diamond to the SH phase is found to 
occur at 44 GPa |Q , which corresponds to an overpres- 
surization of almost a factor of 3. Here we show instead 
that with our new method the transition can be observed 
with affordable computational effort at T = 300 K and 
P = 16 GPa, i.e. very close to the critical pressure. 

A metadynamics simulation was run with the parame- 
ters W = 8.6 eV and Sh = 1 A which are compatible with 
the guidelines given earlier, taking into account that L w 
15 A and a typical Si elastic constant value is c w 100 
GPa. We have preferred a relatively large value for Sh in 
order to enhance volume fluctuations and thus to favour 
the change of volume which accompanies the pressure- 
induced transformation of diamond Si. This choice is 
also instrumental in avoiding that the system makes a 
fake transition to the same crystal structure, since such 
a transition obviously conserves volume. The origin of 
these fake transformations is to be found in the fact that 
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a given crystal structure can correspond to different val- 
ues of h. This problem, which is a consequence of what 
is known as modular invariance |6| , can be fixed in many 
different ways. However, in view of the simple solution 
found here we have not pursued these alternatives, which 
will be discussed elsewhere. 

At each metadynamics step we equilibrated the sys- 
tem for 1 ps and averaged the pressure tensor for an- 
other ps; temperature was controlled by Nose-Hoover 
chains Jl7[. The history of the run is shown in Fig.l. 
The indicator F ■ F g /(|F||F g |) (Fig.l (e)) clearly shows 
that at metastep 35 there is a phase transition (see also 
Fig. 1(d) and Fig.2(a)-(d)). Consequently after this step 
we switched off the Gaussians to let the system evolve 
towards the new structure in a steepest-descent-like man- 
ner till metastep 50 (see Fig.2(e)-(h)). This final struc- 
ture was then evolved for 5 ps with a Parrinello- Rahman 
simulation. During this time very little relaxation of 
the cell parameters took place, which confirms that we 
have reached a minimum of Q(h). A visual inspection 
of the final structure (Fig. 2(h)) as well as an analysis of 
the diffraction peaks showed that the system had made 
a transition to the SH structure, whose parameters are 
a = 2.61 A and c = 2.48 A. We also calculated the atomic 
volume, shown in Fig. 1(c). After step 34 we observe a 
pronounced drop, agreeing well with the results of the 
T = calculation which predicts a change from 17.5 A 3 
per atom in the diamond to 14.8 A 3 per atom in the SH 
phase. 

Another run was carried out in which after metastep 35 
the Gaussians were kept switched on. This metadynam- 
ics led to a series of transitions between the stable SH and 
metastable /3-tin structures with different orientations. 
We also performed a simulation of decompression of the 
SH structure at p = 5 GPa which after 48 metadynamics 
steps transformed into a tetrahedrally-coordinated amor- 
phous structure. Experimentally Si upon decompression 
also does not come back from the /3-tin to the diamond 
structure but transforms to a series of metastable struc- 
tures |lq| . These two examples demonstrate that the 
method is able to find also metastable crystalline and 
amorphous phases. 

This calculation shows that the method overcomes 
many of the limitations of the previous approaches. It 
must be stressed that, in contrast to previous work|Q, ^, 
||, |[ ||, [| 0, |[ ||, [To), it is not a constant-pressure sim- 
ulation method but rather a method for exploring the 
dependence of the Gibbs free energy on the h variables. 
The use of the history-dependent metadynamics allows 
large energy barriers to be overcome in a short time and 
makes this approach very efficient. If other internal slow 
degrees of freedom besides the h variables are present, 



as in molecular crystals, these can naturally be added 
and taken into account within the general scheme [^3). 
Metadynamics is in principle able to visit any crystalline 
structure that is at least metastable at a given pressure; 
however, use of pressures considerably different from the 
critical one may result in the need for a longer time in 
order to escape from the initial minimum. In the exam- 
ple presented here the total aggregated simulation time 
is about 100 ps. This makes the method suitable also 
for an ab-initio MD simulation of systems containing 
about 100 atoms. Another advantage is that ab-initio 
constant-volume codes can be used, avoiding the need to 
use expensive tricks to deal with the Pulay correction |Q. 
In conclusion it can be confidently stated that this new 
method, with its ability to induce structural transitions 
at equilibrium conditions, can substantially improve the 
predictive power of solid-state simulations. 
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FIG. 1: Evolution of cell lengths (a), cell angles (b), atomic 
volume (c), selected peaks of the diamond (averaged over all 
equivalent directions, orange) and SH (blue) structure factor 
(d) and relative orientation of forces F and F 9 (e) during the 
metadynamics. Note the structural transition at step 35. The 
Gaussian term in Eqn.(|^) is switched off after step 35. The 
light blue curve in (e) shows the continuation of the orange 
run in the mode in which the Gaussians are added at every 
metadynamics step. 
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FIG. 2: (a) - (d) Evolution of atomic configurations during 2 
ps of microscopic dynamics (at intervals of 667 fs) across the 
transition at metastep 35. The initial diamond structure (a) 
is strongly strained, compressed along one axis and elongated 
along perpendicular ones (see also Fig.l (a)). In the next two 
snapshots (b), (c) the gradual disappearance of the diamond 
structure can be observed; at the same time, a new periodic 
structure emerges (d). (e) - (h) Evolution during 15 subse- 
quent steps of metadynamics. Note the gradual formation of 
the simple hexagonal phase. From the analysis of diffraction 
peaks the final supercell (h) was found to contain 222 prim- 
itive cells; therefore 6 vacancies are actually present in the 
structure. 



